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Abstract 



We present a recursive procedure to calculate the parameters of the 
recently introduced multicanonical ensemble and explore the approach 
for spin glasses. Temperature dependence of the energy, the entropy 
and other physical quantities are easily calculable and we report re- 
sults for the zero temperature limit. Our data provide evidence that 
the large L increase of the ergodicity time is greatly improved. The 
multicanonical ensemble seems to open new horizons for simulations 
of spin glasses and other systems which have to cope with conflicting 
constraints. 
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The theoretical understanding of spin glasses, for a review see has re- 
mained a great challenge. In particular the low temperature limit leaves many 
open questions about the effects of disorder and frustration. For instance, 
it has remained controversial whether Parisi's [[J mean field theory provides 
the appropriate description for 3D spin glasses. The attractive alternative is 
the droplet model ||, which in turn is equivalent to a one parameter scaling 
picture [[|. The simplest spin glass system to study such questions numeri- 
cally is the Edwards- Anderson model. In its Ising version it is described by 
the Hamiltonian 

H = - JijSiSj, (1) 

<ij> 

where the sum goes over nearest neighbours and the exchange interactions 
Jij = ±1 between the spins Sj = ±1 are quenched random variables. In 
our investigation we impose the constraint Jij — for each realization. 
Recent simulations || of the 3D model in a magnetic field support the mean 
field picture. However, one may argue that sufficiently low temperatures on 
sufficiently large systems have not been reached. For previous simulations of 
the Edwards- Anderson model without magnetic field see || |7| . 

Low temperature simulations of spin glasses suffer from a slowing down 
due to energy barriers. To illustrate the problem, let us consider a simple 
ferromagnet: the 2D Ising model on an 50 x 50 lattice. In Figure 1 we give 
its magnetic probability density versus (3 = T" 1 . The two distinct branches 
below the Curie temperature are associated with free energy valleys in con- 
figuration space, each of which defines a (pure) thermodynamic state. At 
temperatures below the Curie point the ergodicity time 1 r| increases expo- 
nentially fast with lattice size, asymptotically like exp[/ s (/3)L D ~ 1 ], where f s 
is the surface free energy. Therefore, on large lattices at sufficiently low tem- 
perature the simulation of the system will, given a reasonable finite amount of 
computer time, never tunnel from one phase to the other. Besides for partic- 
ular problems (like the determination of the order-order surface free energy) 
this lack of tunneling does not pose a major handicap for Ising model sim- 
ulations. The reason is that the two configuration space valleys are related 
by the exact symmetry Sj — > —Si of the Hamiltonian. Exploring one valley 

is for the present investigation of spin glasses more appropriate to use the term 
ergodicity time r e , instead of tunneling time r* which is appropriate in the context of 
surface tension investigations. 
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Figure 1: Ising model magnetic probability density from 50 x 50 lattice. 
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by means of a simulation yields also all the properties of the other one and, 
hence, allows to overlook the entire system. 

The situation is much more involved for spin glasses. For low enough 
temperature the system is supposed to split of into many thermodynamic 
states, separated by similar tunneling barriers as the two pure states of the 
Ising model. However, unlike in the Ising model the states are not related 
to each other by a symmetry of the Hamiltonian. Rather they appear be- 
cause of accidental degeneracy which in turn occurs because of randomness 
and frustration of the system. For computer simulations this means that one 
would like to explore many independent configuration space valleys while 
keeping track of their relative weights. The groundstate energies associated 
with these valleys may or may not be degenerate, but it should be noted 
that even if they are not degenerate, tunneling between the valleys would 
still be governed by the energy barriers. The physics of these barriers is 
far less well understood as in the ferromagnetic case. As detailed finite size 
scaling (FSS) studies do not exist, it is unclear to us to what extent these 
barriers depend on the system size, whereas the temperature dependence has 
been investigated |]]. We use the notation bifurcation temperature (bifurca- 
tion point) for the temperature at which the spin glass configuration space 
(phase transition or not) begins to split of into a number of valleys which are 
well separated by energy barriers. In the present paper we suggest that the 
increase of r£ can be reduced to a fairly decent power law by performing a 
simulation which covers in a single ensemble a whole temperature range from 
well above to far below the bifurcation point. The appropriate formulation is 
provided by a generalization of the multicanonical ensemble ||. To test the 
approach we have performed multicanonical simulations of the 2D Ising fer- 
romagnet and then of the 2D Edward- Anderson Ising spin glass model. We 
concentrate on ground state properties what is a kind of worst case scenario 
for the performance of the multicanonical simulation. It should be noted that 
Figure 1 does rely on a multicanonical simulation, what explains the slight 
asymmetry between the two branches. 

Ever since the pioneering paper by Metropolis et al.|J most MC simula- 
tions concentrated on importance sampling for the canonical Gibbs ensemble. 



It has always been well-known, for instance |10| , that one is allowed to choose 



phase-space points according to any other probability distribution, if it is con- 
venient. However, a systematic reasoning for a better than canonical choice 
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has rarely been put forward 1 . It is our suggestion that in a large class of 
situations, in particular those where canonical simulations face severe ergod- 
icity problems, it is more efficient to reconstruct the Gibbs ensemble from 
a simulation of a multicanonical ensemble || than simulating it directly. In 
canonical simulations configurations are weighted with the Boltzmann factor 
Pb(E) = exp(—/3E). Here E is the energy of the system under consideration. 
The resulting canonical probability density is 

P C {E) ~ n(E)P B {E), (2) 

where n(E) is the spectral density. In order of increasing severity problems 
with canonical spin glass simulations are: 

i) Simulations at many temperatures are needed to get an overview of the 
system. 

ii) The normalization in equation (2) is lost. It is tedious to calculate im- 
portant physical quantities like the free energy and the entropy. 

Hi) The low temperature ergodicity time r£ diverges fast with lattice size 
(either exponentially or with a high power law). The relative weights of pure 
states can only be estimated for small systems. 

Let us choose an energy range E min < E < E max and define for a given 
function /3(E) the function a(E) by the recursion relation (with the Hamil- 
tonian (1) the energy changes in steps of 4) 

a(E - 4) = a(E) + \fi(E - 4) - (3(E)] E, a(£ max ) = 0. (3) 

The multicanonical ensemble [Q is then defined by weight factors 

P M (E) = exp [-/3(E)E + a(E)}, (4) 

where /3(E) is determined such that for the chosen energy range the resulting 
multicanonical probability density is approximately flat: 

P m u(E) = c mu n(E)P M (E) w const. (5) 

In the present study we take -E max = ((3(E) = for E > E max ) and 
E min = E° the ground state energy of the considered spin glass realization. 

1 Notable may be microcanonical simulations, but for the ergodicity problems on which 
we focus here they perform even worse than the canonical approach does. 
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A multicanonical function /3(E) can be obtained via recursive MC calcu- 
lations. One performs simulations (3 n (E), n = 0,1,2,..., which yield prob- 
ability densities P n (E) with medians -^median- For E < E^ hl < -E™ cdian the 
probability density P n (E) becomes unreliable due to insufficient statistics, 
caused by the exponentially fast fall-off for decreasing E. We start off with 
n = and (3°(E) = 0. The recursion from n to n + 1 reads 



(3 n+1 (E) 



{ (3 n (E)^ for E > Median! 
ft (-^median) 

(££ edian - E)- 1 In [P n (E^ eAi J/P n (E)\ for ££ edian > £ > ££ ir 
I P n+ \E- in ) for £ < ££ in . 

(6) 

Here the n th simulation may be constrained to E < -^median by rejecting 
all proposal with energy E > ^median > but one has to be careful with such 
bounds in order to maintain ergodicity. The recursion is stopped for m 
with E™^ 1 = E° being groundstate. Starting with this simple approach we 
have explored several more sophisticated variants. Considerable speed-ups 
and gains in stability could be achieved. The CPU time spent to estimate 
the multicanonical parameters was 10% to 30% of the CPU time spent for 
simulations with the final set. 

Once the functions /3(E) and a(E) are fixed, the multicanonical simula- 
tion exhibits a number of desireable features: 

i) By reweighting with exp[— (3E + f3(E)E — a(E)] the canonical expectation 
values 

00) = Z0)~ 1 J2O(E) n(E) ex V (-(3E), (7) 

E 

where Z(/3) = Y,E n (E) exp(—/3E) is the partition function, can be recon- 
structed for all (3 in an entire range /9 min < (3 < /3 max . Here (3 min = (3(E max ) 
and /? max = (3(E min ) follow from the requirement -E max > E((3) > E min , and 
E(/3) follows from (7) with 0(E) = E. With our choice E max = and 
E m[n = E° groundstate, /5 m i n = and /3 max = oo follows. 

ii) The normalization constant c mu in equation (5) follows from Z(0) = 
Y;E n (E) = 2 N , where N is the total number of spin variables. This gives 
the spectral density and allows to calculate the free energy as well as the 
entropy. 

Hi) We conjecture that the slowing down of canonical low temperature spin 
glass simulations becomes greatly reduced. For the multicanonical ensemble 
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it can be argued [§] that single spin updates cause a ID random walk be- 
haviour of the energy E. As E max — E min ~ V, one needs V 2 updating steps to 
cover the entire ensemble. For first order phase transition the observed slow- 
ing down was only slightly worse than this optimal behaviour. Our present 
MC data show more drastic modifications for spin glass simulations. 

To quantify our discussion of the slowing down, we have to define the 
ergodicity time r£. Roughly speaking it is the CPU time needed to collect 
independent configurations of the system, for instance groundstates when 
this is the main interest. Regarding the definition of r£, an L-independent 
over-all factor is free. We define r£ as the average number of sweeps needed 
to move the energy from E max to E min and back. A sweep is defined by up- 
dating each spin on the lattice once (in the average). It should be noted that 
the number of "tunneling" events with respect to the ergodicity time gives 
a lower bound on the number of independent groundstates samples. This 
follows from another trivial, but remarkable property of the multicanonical 
ensemble: Each time a sweep is spent at /3(E) = 0, the memory of the pre- 
vious Markov chain is lost entirely, and a truly independent new series of 
configurations follows. The condition E > E max is appropriate to substitute 
for the somewhat too strict constraint of a sweep at (3 = 0. As a corol- 
lary: with a disordered starting configuration the multicanonical ensemble is 
immediately in equilibrium. 

As an exercise and to check our code on exact results, we performed a 
multicanonical simulation of the 2D Ising model with < (3 < oo. We 
kept the time series of two million sweeps and measurements on a 25 x 25 
lattice and verified that the finite lattice specific heat results of Ferdinand 
and Fisher |12| are well reproduced. No difficulties are encountered with the 
multicanonical ensemble when crossing the phase transition point. To explore 
the possibility of zero temperature entropy calculations, we used Z(0) = 2 625 
as input and obtained S° = 0.61 ± 0.09 for the total groundstate entropy. 
This corresponds to an estimate of 1.84 ± 0.17 groundstates, i.e. within 
statistical errors we are in agreement with two. Using Z(0) = 2 2500 and a 
time series of four million sweeps on a 50 x 50 lattice we obtained 2.07 ±0.22 
for the number of groundstates. Figure 1 is obtained from the simulation of 
this lattice. 

After this test we turned to the 2D Edwards-Anderson spin glass. On 
lattices of size L = 4, 12, 24 and 48 we performed multicanonical simula- 
tions. Up to L = 24 we investigated ten different realizations per lattice and, 
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Figure 2: Multicanonical energy density distribution for one L = 48 realiza- 
tion. 



due to CPU time constraints, we considered only five realizations for the 
L = 48 lattice. The multicanonical energy distribution for one of our L = 48 
realizations is depicted in Figure 2. The fall-off for — e < is like that of 
the canonical distribution at $ — 0. For < — e < — e° an impressive flat- 
ness (about 800 energy entries on the lattice under consideration) is quickly 
achieved by the recursion (6). Close to the groundstate some fluctuations are 
encountered on which we comment elsewhere [13]. As it is not obvious from 
the scale of the figure, we like to remark that the groundstate is not the state 
with the lowest number of entries, but a state close to it. Further, it should 
be understood that deviations from the desired constant behaviour (5) do 
only influence the statistical error bars, but not the estimates themselves. 
Therefore, such deviations do not pose problems as long as they can be kept 
within reasonable limits of approximately one order of magnitude. 
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Figure 3: Ergodicity times versus lattice size on a double log scale. 
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The Table gives an overview of some of our numerical results. The er- 
godicity time r£ is as defined above, and for /3 max we take /3(E°), where it 
should be noted that due to our computational procedure /3(E) is a noisy 
function. The reported values and their error bars are obtained by combining 
the results from the different realizations, which enter with equal weights. In 
Figure 3 we plot the ergodicity time versus lattice size L on a log-log scale. 
The data are consistent with a straight line fit (Q denotes the goodness of 
fit), which gives the finite size behaviour 

T e _ £4.4(3) sweeps . (g) 

In CPU time this corresponds to a slowing down ~ y 3 2 ( 2 ). It should be 
remarked that a fit of form r£ ~ exp(cL) results in a completely unacceptable 
goodness of fit Q < 1CT 6 . Still, the behaviour (8) is by an extra volume factor 
worse than the close to optimal performance we hoped for. The reasons will 
be considered elsewhere [ fT3"| . 



L 


Statistics 


rt 


/^max 


4 


10 x 1,600,000 


35.3 ±2.8 


0.75 ±0.06 


12 


10 x 760,000* 


2, 607 ±450 


1.47 ±0.04 


24 


10 x 1,600,000 


193, 750 ±43, 820 


2.12 ±0.11 


48 


5 x 5,200,000* 


1,457, 315 ±516, 925 


2.22 ±0.07 



Table 1: Overview of some results. For the data points marked with * the 
statistics for different realizations varies somewhat and average values are 
given. 



We estimate the infinite volume groundstate energy and entropy from 
FSS fits of the form fi = foo + c/V. The entropy fit is depicted in Figure 4, 
and the energy fit looks similar. Our energy estimate is e° = —1.394 ±0.007, 
consistent with the previous MC estimate M e° = —1.407 ± 0.008 as well 



as with the transfer matrix result [14] e = —1.4024 ± 0.0012. Our entropy 



estimate s° = 0.081 ± 0.004 is also consistent with the MC estimate [|7J 
s° = 0.071 ± 0.007, but barely consistent with the more accurate transfer 
matrix result |14j s° = 0.0701 ± 0.005. Still, a larger statistical sample 
and a careful study of systematic error sources would be needed to claim 
that there is a significant discrepancy. Our results could be improved by 
exploiting the high temperature expansion of the entropy, as it was done in 
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Figure 4: FSS estimate of the infinite volume entropy per spin. 
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M. However, this would be against the spirit of this paper which is to explore 
the possibilities and limits of multicanonical spin glass simulations. 

It is not entirely straightforward to compare multicanonical and standard 
simulations. For instance autocorrelation times of multicanonical simulations 
come out short due to the triviality that the simulation spends most of its 
time at rather small effective (3 values. Therefore, autocorrelations are not 
well suited. Our ergodicity time, the average number of sweeps to find truly 
independent groundstates, is a more useful quantity. When relying on it, 
i.e concentrating on groundstate properties, we should remember that this 
pushes the multicanonical simulation to its extreme. This way of calculating 
groundstate properties gives for free all properties in-between, from f3 = 
on. If a phase transition occurs, its study is also included (we noticed no 
particular difficulties in the 2D Ising model). For instance in the 3D Ising 
spin glass the canonical slowing down at (3 C JTJ is already worse than the one 
given by our equation (8). 

Although the slowing down (8) is severe, it seems to provide an impor- 
tant improvement when compared with the slowing down which canonical 
simulations encounter for temperatures below the bifurcation temperature. 
For L > 24 canonical simulations [[[], Q are unable to equilibrate the sys- 
tems at the /3 max values reported in our table since the relaxation time 
is by far too long. Presumably due to this fact the literature focused on 
other questions and does not provide detailed FSS investigation of relax- 
ation times. At the present stage we did not want to spent CPU time on 
extensive canonical simulations. Therefore, detailed comparisons are not 
possible. However, a rough estimate of the canonical ergodicity time may 
be i r c e anonical ~ exp (c$ - C') with C « 15 and C » 11.6. The 
scaling with L may then be hidden in the L dependence of our /3 max values, 
which is argued to be divergent like /3 max ~ hi(V), and this line of reasoning 
gives a slowing down of the canonical algorithm like V 15 . With /9 max = 2.12 
(our L = 24 case) one gets a canonical ergodicity time of order 10 8 — 10 9 
when the missing constant is assumed to be of order one, whereas our r| 4 is 
about 2 x 10 5 . Already for moderately large /3-values, a doubtful procedure 
has emerged in the literature [§, [5]: Realizations which (according to a well 
defined criterium) cannot be equilibrated are simply omitted. As these real- 
izations £1X6 cLS legitimate members of the statistical ensemble as any other 
choice of the quenched random variables Jij, an uncontrollable bias is in- 
evitably introduced. In contrast to this, we were able to equilibrate every 
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single realization encountered. 

When one is only interested in groundstate properties, minimization algo- 
rithms have to be considered. As a method simulated annealing [IT| stands 
out because of its generality, although there are more efficient algorithms for 
special cases, which should be used when appropriate. In simulated annealing 
the results depend on the cooling rate r = — AT/sweeps. For our model the 
behaviour e(r) = e° + c r 2 with c ~ 0.5 (AT = —0.1 fixed) is indicated [|15|1 . 
To find a true groundstate, one has to reduce [e(r) — e ] to the order l/V. 
Assuming that the constant c is volume independent (only the lattice size 
100 x 100 was considered in [|HJ), this translates to # sweeps ~ V 4 = L 8 , 
far worse than our equation (8). This result is kind of amazing as the multi- 
canonical ensemble has eliminated directed cooling and is nevertheless more 
efficient. If one does not insist on true groundstates one can then relax the 
condition [e(r) — e°] ~ l/V". For instance, any behaviour [e(r) — e°] — >• with 
L —>■ 00 would still give the correct density, and simulated annealing would 
slow down less dramatic. On the other hand, this would also imply a less 
stringent multicanonical simulation. 

A comparison with the cluster-replica MC algorithm |7| is even less clear 
cut. The obtained estimates of the groundstate energy and entropy are in 
accuracy similar to ours. As one has to simulate many replica at many (3- 
values a direct comparison is impossible. Clearly, the results reported on 
slowing down are more promising than ours for the ergodicity time. On 
the other hand, to our knowledge the replica MC approach has never been 
applied to the 3D Edwards-Anderson model, and one may well encounter 
difficulties. In contrast, 3D multicanonical simulations are straightforward. 
In fact the dimension is just a parameter in all our computer programs and 
we have already carried out various 3D test runs. Let us address the ques- 
tion of algorithms from a more general perspective. With a Metropolis type 
implementation our optimum performance will be bounded from below by 
a slowing down ~ V 2 in CPU time, and we are outperformed by any algo- 
rithm which can do better than this. For a number of important, but often 
highly specialized, applications such better algorithms exist and should be 
used. The main advantage of the multicanonical ensemble is its general- 
ity. With this respect our method resembles simulated annealing JTT|, while 
clearly avoiding some of its disadvantages. Most notably, the relationship 
to the equilibrium canonical ensemble remains exactly controlled. Finally, it 
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should be stressed that the multicanonical ensemble is an ensemble and not 
an algorithm. One may find better algorithms than conventional Metropolis 
updating to simulate this new ensemble. To try a combination with cluster 
algorithms is an attractive idea. 

Our results make clear that the multicanonical approach is certainly a rel- 
evant enrichment of the options one has with respect to spin glass simulations. 
The similarities of spin glasses to other problems with conflicting constraints 



1 1] suggest that multicanonical simulations may be of value for a wide range 



of investigations: optimization problems like the travelling salesman, neu- 
ral networks, protein folding, and others. Multicanonical simulations of the 
3D Edwards- Anderson model may eventually shed new light on the question 
whether the model exhibits mean field like behaviour or some kind of droplet 
picture applies p], [|, f|, [5[]. It seems that previous numerical work remained 
somewhat inconclusive as sufficiently low temperatures could not be reached 
without destroying the thermodynamic equilibrium. 
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